Radial distribution function of semiflexible polymers 



\0 . 
0^ ■ 

<^ : 

Oh! 

<D . 
00 ■ 

l> ■ 



> 

o 
in 

o 

ON 



o 



Jan Wilhelm and Erwin Frey 
Institut fiir Theoretische Physik, Technische Universitat Miinchen, 85747 Garching, Germany 

(Phys. Rev. Lett. 77, 2581 (1996)) 

We calculate the distribution function of the end-to-end distance of a semiflexible polymer with 
large bending rigidity. This quantity is directly observable in experiments on single semiflexible 
polymers (e.g., DNA, actin) and relevant to their interpretation. It is also an important start- 
ing point for analyzing the behavior of more complex systems such as networks and solutions of 
semiflexible polymers. To estimate the validity of the obtained analytical expressions, we also deter- 
mine the distribution function numerically using Monte Carlo simulation and flnd good quantitative 
agreement. 



While we have a comparably complete theoretical pic- 
ture of highly flexible chain molecules, the statistical me- 
chanics of semiflexible polymers is a field with a number 
of open questions that has received renewed attention 
lately. Considerable motivation stems from the crucial 
importance of the elasticity of biopolymers like spectrin, 
actin, and microtubules for the mechanical properties of 
cells Q. Recent advances in visualizing and manipu- 
lating such macromolecules have provided unique exper- 
imental tools for the study of the static and dynamic 
properties of single filaments [^-||. A quantitative mea- 
sure for a polymer's flexibility is its persistence length £p, 
which is the characteristic length governing the decay of 
tangent -tangent correlations. Nature provides polymers 
of very different stiffness, e.g., £p « 10 nm for spectrin 0, 
ip K. 50 nm for DNA §], ^j, « 17/Ltm for actin [|,§ and 
Ip « 5.2 mm for microtubules On length scales larger 
than a few £p a polymer with contour length L 3> €p 
(flexible polymer) can be described as a self-avoiding 
freely jointed chain. For molecules with t := L/lp of 
the order of one (semiflexible polymer) this is not possi- 
ble, as can be seen immediately by comparing a typical 
contour (e.g., from Ref. |^) to the random walk corre- 
sponding to a freely jointed chain. Even for long strands 
of DNA with t K, 1000 measurements of the molecule's 
extension at large forces have revealed significant devia- 
tions from the freely jointed chain model which find 
their explanation in bending elasticity . Thus it is es- 
sential to consider models which take chain rigidity into 
account. For sufficiently stiff chains effects resulting from 
self-avoidance can be neglected due to the strong ener- 
getic suppression of configurations where the chain folds 
back onto itself. Furthermore the polymers under con- 
sideration can be regarded as inextensible The cor- 
responding model is the wormlike chain introduced by 
Kratky and Porod almost 50 years ago pO[ |. 

A central quantity for characterizing the conforma- 
tions of single polymer chains is the distribution func- 
tion G(r; L) of the end-to-end distance r for given con- 
tour length L and persistence length ip. For models like 
the wormlike chain with only short-range interactions 
between monomers, G(r; L) actually is the probability 
density of finding any two monomers at relative position 



r = r(s) — r(s') where L — \s — s'\ is the distance be- 
tween the monomers along the chain. For a freely jointed 
phantom chain G(r; L) is known exactly Jill and for 
many purposes is well approximated by a simple Gaus- 
sian. Complications arise when the self-avoidance of real 
chains is taken into account |12|. In previous investi- 
gations of the wormlike chain G(r; L) was obtained ap- 
proximately for almost fully flexible polymers (large t) in 
the form of corrections to the Gaussian distribution func- 
tion up to order t^"^ [ p^p4| . For general t only the low- 
est three even moments were calculated analytically . 
Higher even moments were obtained by numerical tech- 
niques . For polymers close to the stiff limit all even 
moments of the distribution function were calculated in 
an expansion in t |l5| , |l7| , but the expressions obtained 
for G(r; L) in this limit are only given up to quadratures 
|l7| , ^ and do not show the correct qualitative behavior 
when integrated numerically. 

In this letter we determine G(r; L) in two- and three- 
dimensional embedding space in an approximation valid 
for small t and compare the analytical expressions ob- 
tained to data from a Monte Carlo simulation. The range 
of validity of our results almost extends to values of the 
bending rigidity where the Daniels approximation 14 
becomes applicable. 

For our analytical calculations we adopt a continuum 
version of the wormlike chain model where the polymer is 
represented by a differentiable space curve r(s) of length 
L parameterized to arc length |l^. Its statistical prop- 
erties are determined by an effective free energy 
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where t(s) = dr{s)/ds is the tangent vector at arc length 
s. The resulting persistence lengths are £p = n/ksT 
for d = 3 and ip = 2K/kBT for d = 2, where d is the 
dimension of the embedding space [ pO[ . The inextensi- 
bility of the chain is expressed by the local constraint 
|t(s)| = 1 which leads to non-Gaussian path integrals. 
Note that this local curvature model is equivalent to a 
one-dimensional nonlinear cr-model. While a few quanti- 
ties like (i?^) and (i?^) can be obtained exactly, one has 
to resort to some approximative scheme to calculate the 
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end-to-end distribution function 



G(r;L) = (5(r-R)) 



(2) 



where R := r{L) — r(0) and 6{r) is the Dirac (5-function. 

In order to understand the effect of possible approxi- 
mations to Eq. (H) for stiff polymers, it is instructive to 
recapitulate the classical problem of bending a rigid rod. 
The energy of a straight rod of length L and bending 
modulus K is an almost linear function of its end-to-end 
distance r: Ec\ ~ fc{L — r) where fc = ktt^/L^ is the crit- 
ical force for the onset of the Euler instability (since the 
radial distribution function docs not specify the direction 
of the tangent vectors at s and s — L, the appro- 
priate boundary conditions are open ends). Neglecting 
fluctuations around the classical contour this would lead 
to an end-to-end distribution function with maximum 
weight at r = L, G(r; L) oc exp[-/c(L - r)/{kBT)]. But, 
for a completely stretched chain there is up to global ro- 
tations only one possible configuration and consequently 
the end-to-end distribution function has to vanish at full 
extension. Hence it is essential to take into account en- 
tropy effects. While stiff chains are energy dominated, 
more flexible chains are mostly governed by the entropic 
effects with the limit being the freely jointed chain. Ap- 
proaching full extension, however, G'(r; L) must vanish 
for all chains. This shows that while G{r; L) tends to 
the distribution 5{L — r)/47rL^ [d = 3) for k oo and 
t — !■ 0, it can never be expanded in t around this limit 
(consider the probability of full extension). For this rea- 
son, it is impossible to obtain the distribution function 
from an expansion of all even moments in terms of t . 
Softening the constraint of fixed contour length affects 
both energetic and entropic contributions to the distri- 
bution function in essential ways. If it is relaxed to the 
point of fixing only {^^ (ist^(s]) by means of a single La- 
grange multiplier (e.g., Ref. the distribution func- 
tions obtained are essentially Gaussian and will not show 
the correct qualitative behavior for stiff polymers. Fail- 
ure to reproduce the vanishing of G(r; L) near full ex- 
tension also results when approximations are used which 
neglect essential parts of the fluctuational contributions 
as in Ref. @. 

For end-to-end distances r close enough to full stretch- 
ing and/or large values of k, the typical configuration 
of the chain will be close to a straight rod. Thus the 
deviations of the tangent vectors from the average di- 
rection can be treated as small variables. For c? = 3 
we parameterize the contour through the tangent field: 

t(s) = {ax{s),ay{s), 1)/^^! + 0^(3) -I- a2(s), which prop- 
erly takes into account the constraint of inextensibility. 
We employ a harmonic approximation and keep only 
terms up to second order in 7i, the measure factor, and 
the arguments of the (5-function in Eq. (|2|). The error 
caused by this approximation vanishes near full exten- 
sion. From here on we shall measure all lengths in units of 
L and all energies in units of ksT . With this convention 
we have (for c? = 3) k = £p and t = k"-*^ . We also drop the 
second argument in G(r; 1) and make rotation invariance 



expHcit by writing G(r). Use of 2it6{x) = J dqexp{iqx), 
expansion of ax{s) and ay{s) in cosine series as appropri- 
ate for the boundary conditions of open ends, and eval- 
uation of the resulting Gaussian integrals leads to 



G(r) 



2k 



47r7V^ 



(3) 



where A/" is a normalization factor compensating the fail- 
ure of the approximation used to conserve the normal- 
ization of G(r) for finite k. Details of the calculation 
can be found in a forthcoming publication [ p2[ . For 
k(1 — r) ^ 0.2 the distribution function is dominated 
by the fc = 1 term which is just our heuristic result from 
above ( "Euler instability" ) if the dimensional factors are 
put back in. For r 1, however, more and more terms 
have to be considered. By writing Eq. (^ as an inte- 
gral over Fourier expanded ^-functions it is possible to 
transform Eq. (0) to 
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where H2{x) = Ax'^ — 2 is the second Hcrmite polyno- 
mial. This series converges very quickly for k(1 — r) ^ 0.2 
where the behavior of G(r) is completely dominated by 
the £ — 1 term. 

In recent experiments as well as in simulations the 
polymer was effectively restricted to a d = 2 embed- 
ding space ||^-|6|,p3|. In this case it is convenient to 
choose a parameterization t(s) = (cos sin 0). The 
resulting effective free energy is quadratic in 0, 7i = 
(k/2) / ds {d(l)/dsf . In order to calculate G(r) we again 
approximate the constraints given by the i5-functions 
in Eq. (H) by = ry{l) « / dscf) and r = rx{l) « 
1 — i J ds4)'^ . A calculation along similar lines as before 
[E2i leads to the result 
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with D-^i2{x) a parabolic cylinder function. The conver- 
gence properties of Eq. (||) are similar to those of Eq. (Q). 

In order to assess the quality of our approximations, we 
have used Monte Carlo simulation to evaluate G(r) nu- 
merically. We adopted the following discretized version of 
the wormlike chain: The polymer is described as a chain 
composed of N tethers of fixed length a — 1 /N and direc- 
tion t with a bending energy TLj^ = ka~^ Si=i^ ' ^i+i 
p^ . The standard Metropolis algorithm was used to 
measure G(r). We found that results cease to depend 
appreciably on N as soon as there are three to four seg- 
ments in one persistence length. On the order of 10^ MC- 
steps per segment were performed. Final results were ob- 
tained by averaging over several independent runs. The 
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accuracy of (R) obtained was typically of the order of 
0.5%. Measured expectation values (i?^) and (i?"*) were 
in agreement with known exact expressions up to the es- 
timated statistical errors. 
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FIG. 1. Comparison of G{r) from the Monte Carlo simu- 
lation (symbols) for t = 10, 5, 2, 1 and 0.5 (left to right) to 
Eq. (^) (solid lines). The dashed and dot-dashed lines show 
the second Daniels distribution for t = 10 and t = 2. For 
t = 10 it almost coincides with the numerical data while it is 
quite far off for t — 2. Error bars for the Monte Carlo data 
are approximately of the size of the symbols. 



Fig. shows a comparison of the normalized G(r) from 
Eq. (|^) to the data from our Monte Carlo simulation. 
Note that there is no free parameter to adjust the curves. 
The curves for the d = 2 case are qualitatively similar and 
the agreement between theory and the d = 2 Monte Carlo 
simulation is of equal quality. It can be improved some- 
what further by applying a simple correction procedure 
specific to the d 
is that Eqs 

lively right and are good approximations for ^.j. 




221. The general observation 
) reproduce the data qualita- 
ip^ 0.5i, 

that is t ^ 2. However, even for large values of k one 
would not expect the harmonic approximation to yield 
acceptable results for small r since the curvature of the 
polymer must then be large. The somewhat surprising 
quality of the approximation in this region can be at- 
tributed to the fact that there the distribution function 
is dominated by the energy Ec\{r) of the most probable 
configuration. As discussed below Eq. (^, the dominant 
linear term of Ec\{r) is reproduced by the harmonic ap- 
proximation, which corresponds to the well known fact 
that /c can be obtained by looking only at infinitesimal 
compressions. For r « 0, however, the deviations of E^^i 
from the linear form get significant for all n. While this is 
irrelevant for most applications of the distribution func- 
tion due to the small probability of such configurations, 
it implies that the ring-closure probability G(0) is not re- 
produced correctly by our approach. Instead one should 
expand the configurations around the contour of minimal 
energy as in Ref. p5[ where the resulting expressions for 
G(0) are evaluated numerically yielding results in agree- 
ment with our Monte Carlo data. 

Another possibility to check the validity of the ob- 
tained distribution functions is to compare their moments 
to known results. The moments (i?") can either be calcu- 



lated from the given expressions for G{r) or by applying 
the harmonic approximation directly to the generating 
function (e-^^). The latter method yields an expansion of 
(i?") in t which is seen to be correct to 0{t^) when com- 
pared to the results of Ref. . Calculating the moments 
directly from the i —\ term of Eq. (^) for k(1 — r) < 0.2 
and the k = 1 term of Eq. (||) for /t(l — r) > 0.2 pro- 
duces analytical expressions which can be expanded in 
terms of t only up to a small correction vanishing like 
exp[— const. k]. This leads to expansion coefficients that 
are not rational numbers but differ very slightly from 
the results of the generating function method (by about 
0.1% for t = 1). Calculation from the distribution func- 
tion has the advantage that averages (i?") with arbitrary 
a > —d + 1 can be obtained. The results indicate that 
the expansion derived in |l7|| for a — 2n with n = 0, 1, . . . 
is valid for all a > — d -|- 1 at least up to 0{t^). 

A quantity of immediate experimental interest is the 
force-extension relation. For sufficiently stiff polymers it 
can be obtained from the given G{r) for arbitrary forces 
in situations where the ends of the polymer can rotate 
freely. The well known strong- force limit (e.g., |^) is 
reproduced by our distribution functions for arbitrary 
stiffness since it is determined by the behavior of G{r) 
for r — > 1 where the harmonic approximation is expected 
to be good for all values of k: For large forces / we have 
r « 1 and consequently G(r) will be dominated by the 
£ = I term of Eq. (^). The relevant free energy is thus 
F = - log G(r) + /(I - r) « 1/(4k(1 - r)) + /(I - r), 
where we neglected logarithmic terms since 1 — r ^ 
for / — > oo. Since the corresponding distribution func- 
tion will be strongly peaked for large /, (r) is given by 
the position of the minimum of F: (r) = 1 — I/^/AkJ 
which agrees with Ref. |^ . Force-extension relations for 
small forces can be obtained from moments of the dis- 
tribution function. Our results are thus consistent with 
existing linear response treatments [ p6|j2^ to 0{t^). If 
the polymer's orientation at one end is fixed, the charac- 
ter of the response depends on the direction of the applied 
force and G(r|uo; L) is needed |2^ (notation of Ref. |p8|). 
Note that the linear response can be evaluated exactly for 
arbitrary values of t in this case . 

Using fluorescence microscopy, it is possible to visual- 
ize the thermal undulations of single actin filaments con- 
strained to quasi two-dimensional configurations 
Both contour lengths and distances in embedding space 
can be measured from the resulting images. If locality 
of the interactions along the chain is assumed, differ- 
ent segments of one physical polymer of length Lq are 
statistically independent. One can thus probe different 
lengthscales by obtaining experimental distribution func- 
tions G{r; L) for different L < Lq. Using only £p as free 
parameter, our theoretical results can be fltted to the 
experimental distributions. The quality of the fits for 
different L will indicate at what scales actin can actually 
be described by the wormlike chain model. This might 
reveal new interesting physics and help to clarify some 
of the ambiguities that arise in normal- mode analyses of 
actin flickering: The average squared amplitudes of the 
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normal modes with mode number k of the elastic Hamil- 
tonian fail to decay like 1 /fc^ as predicted by the wormlike 
chain model |p|,|6|j23|. Since G{r;L) can be determined 
without differentiation and Fourier transformation of the 
observed contours, the radial distribution function anal- 
ysis is expected to be much more stable against unavoid- 
able experimental errors. Another method that has been 
used to analyze actin flickering is the comparison of the 
measured decay of the tangent-tangent correlation func- 
tion with the exact result (t(s)t(s')) c>c exp(— |s — s'\/£p) 
However, since an exponential decay of correla- 
tions is expected for all systems with only short range 
interactions, this kind of analysis is probably not very 
sensitive to deviations from the wormlike chain model. 
A very interesting possibility would be to attach two or 
more markers (e.g., small fluorescent beads) permanently 
to single strands of polymers and to observe the distri- 
bution function of the markers' separation. This would 
eliminate all the experimental difficulties associated with 
the determination of the polymer's contour. Note that 
contrary to existing methods of analysis it is not neces- 
sary to know the length of polymer between two markers; 
this quantity can be extracted from the observed distri- 
bution functions along with £p. 

In conclusion, we have derived approximate analyti- 
cal expressions for the end-to-end distribution function 
G(r; L) of semiflexible polymers. The essential ingredient 
in our calculation is the choice of a parameterization for 
the polymer's configuration that satisfies the constraint 
of fixed contour length by construction. Comparison with 
Monte Carlo data shows good quantitative agreement 
for t = L/£p ^ 2. Since the Daniels approximation is 
valid for t ^ 10, G{r; L) of the wormlike chain is up to a 
crossover region now available for the whole range from 
rigid rods to random coils. The range of stiffness accessi- 
ble to the approximations used is highly relevant for the 
physics of rather rigid polymers like actin. Knowledge 
of G(r; L) provides new possibilities for the experimental 
determination of tp as well as the lengthscales at which 
real semiflexible polymers can be described by the worm- 
like chain model. Since the actual form of the single chain 
probability distribution function is an important input 
for the theory of many-chain systems, we hope that our 
work will contribute to the investigation of more com- 
plex questions such as the dynamics and viscoelasticity 
of networks and solutions of semiflexible polymers. 
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